%% simulation setting
% sin
bSineQuantization           = 1;        % enable sine signal quantization
    nSineQuantizationBits       = 16;   
    bSineDither                 = 0;    
bTwoShot            = 0;                % add another sinu with Fin2 = 1.5Fin

% pwm
bPwm                = 1;                % enable pwm process
bPwmInterplation    = 0;                % interpolate pwm %-todo
PwmEdgeType         = 0;                % 0: Tail Edge; 1: Lead Edage; 2: Double Edge
PwmLevelType        = 0;                % 0 2-level; 1: 3-level
bPwmPlot            = 0;                % enable pwm plotting



% analyse
bNewFigure          = 0;                % paint in new figures
bHoldOn             = 1;                % hold on and plot with previous curves

bAnalyseSin         = 0;                % enable sin fft analyse
bAnalysePwm         = 1;                % enable pwm fft analyse

bAVG                = 0;
    Navg            = 1;                % 0: no averaging; n: number of averaging

bPwmFftWindowing    = 0;                % windowing for pwm fft
bPwmRemoveDc        = 0;                % remove pwm dc

% snr
Fmax                = 24000;            % max baseband frequency
Fmin                = 100;              % min baseband frequency



%% Constants
A                   = 1;
Ft                  = 6400;                 %*target input frequency, for generation of Fin
Fs                  = 48000;                % input sampling frequency
OSRintp             = 8;                    %*interplator oversampling ratio
Fov                 = Fs * OSRintp;         % oversampling frequency
RESpwm              = 16;                   % pwm resoluation
OSRpwm              = 2 ^ RESpwm;           % pwm oversampling ratio
Fpwm                = Fov * OSRpwm;          % pwm output frequency
Ncycles             = 1;

%% FFT Length Calculation Strategy
%  In case that the calculation complexity exceeds the system's capacity, we need to carefully choose
%  the proper Fin(Ft), Ncycles, Lin, RESpwm and NFFTpwm. Here is intelligent strategy to calculate the
%  values.
Ctotal          = 25;                       % total calculation capacity, by power of 2
Cpwm            = Ctotal;                   % allowed pwm complexity
Cin             = Cpwm - RESpwm;            % allowed input complexity
Lin_c           = 2 ^ Cin;                  % allowed input length
Lcycle          = Fov/Ft;                    % one cycle length of input signal
Ncycles_c       = Lin_c / Lcycle;           % allowed number of cycles
Ncycles         = max(max(primes(Ncycles_c)), Ncycles);   
                                            % real number of cycles
Lin             = Lcycle * Ncycles;         % real input length
Lpwm            = Lin * OSRpwm;             % pwm output length
NFFTin          = 2^nextpow2(Lin);          % NFFT for input signal
NFFTpwm         = 2^nextpow2(Lpwm);         % NFFT for output signal
Fin             = Fov * Ncycles / NFFTin;    % adjust input frequency

%% - no modification by users -
hFig_sin            = 0;
hFig_sin_fft        = 1;
hFig_pwm            = 2;
hFig_pwm_fft        = 3;
% - ------------------------ -

%% initialization
% generated constants
Tov          = 1/Fov;
Tpwm         = 1/Fpwm;
dbref       = -20 * log10(NFFTin/2);


% coherent sampling
%b           = primes(floor(Fin_obj*Nfft/Fs));   % corresponding fft index bin for sinusoid
%Fin         = b(length(b))*Fs/Nfft;             % adjust input frequency for bin
%fprintf(1,'Adjust Input frequency from %2.2f Hz to %2.2f Hz\n',Fin_obj, Fin);


%% generate input signal
t = (0: NFFTin-1) * Tov;
if(bSineQuantization == 0)                      % no quantization
    y = A * sin(2*pi*Fin*t);
    
    if(bTwoShot)
        y2 = A * sin(2*pi*1.5*Fin*t);
        y = (y + y2)/2;
    end
else
    y = sine(nSineQuantizationBits, NFFTin, Fov, Fin, A, bSineDither);          % with quantization
    
    if(bTwoShot)
        y2 = sine(nSineQuantizationBits, NFFTin, Fov, Fin*1.5, A, 0);
        y = (y + y2)/2;
    end
        
end;
%plot(t,y);

%% PWM
tpwm = (0: NFFTpwm-1) * Tpwm;
if(bPwm)
%    [pwm, st, yin] = upwm_v2(y, RESpwm, PwmEdgeType, PwmLevelType, bInterp);
    pwm = upwm_v2(y, RESpwm, PwmEdgeType, PwmLevelType, bPwmInterplation);
end

if(bPwm && bPwmPlot)
    if(bNewFigure) figure; else figure(hFig_pwm); end
    stairs(tpwm,pwm);
  %  hold on;stairs(tpwm, st);
  %  hold on;stairs(tpwm, yin);
end


%% FFT - sin
if(bAnalyseSin)
%   FFT
    x_f_sin = (0:NFFTin-1) * Fov/NFFTin;
    if(nFftAveraging == 0)                         % no averaging
        f_sin = fft(y, NFFTin);
    else                                            % fft averaging with psd2

        f_sin = psd2(y, NFFTin, nFftAveraging);
    end;
    
%   Plot
    sFigName = sprintf('%dBit(s) Sine, Quantization=%d', nSineQuantizationBits, bSineQuantization);
    if(bNewFigure)
        figure('Name',sFigName); 
    else
        figure(hFig_sin_fft); 
        set(hFig_sin_fft, 'Name', sFigName);
    end;
    
    semilogx(x_f_sin, 20 * log10(abs(f_sin)));
    title(sFigName);
    grid on;
end



%% FFT - PWM
if(bPwm && bAnalysePwm)
%   FFT
    x_f_pwm = (0:NFFTpwm-1) * Fpwm/NFFTpwm;
    if(bAVG)
        f_pwm = psd2(pwm, NFFTpwm, Navg);
    else
        if(bPwmFftWindowing)    pwm     = pwm .* blackmanharris(NFFTpwm)';end
        if(bPwmRemoveDc)        pwm     = pwm - mean(pwm); end
        f_pwm = fft(pwm, NFFTpwm);
    end


% SNR
    snr_pwm = snr(f_pwm, NFFTpwm, Fpwm, 3, Fmin, Fmax);
    fprintf(1, 'SNR = %f.2(dB)\n', snr_pwm);
    
% ploting
    sFigName = sprintf('%dBit(s) %d-Level PWM, %d FFT Spectrum, %d OSR, UPWM=%d, SNR=%.2fdB [%dHz-%dkHz]', RESpwm, PwmLevelType+2, NFFTpwm, OSRpwm, bSineQuantization, snr_pwm, Fmin, Fmax/1000);
    if(bNewFigure) 
        figure('Name', sFigName);
    else
        figure(hFig_pwm_fft);
    end
    
    if(bHoldOn) hold on; end
    semilogx(x_f_pwm(1:NFFTin), 20 * log10(abs(f_pwm(1:NFFTin))));
    
    %hold on;plot([Fmin Fmin], [0 120], 'Color', 'red');
    %hold on;plot([Fmax Fmax], [0 120], 'Color', 'red');
    
    title(sFigName);
    grid on;

end;


